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1. Introduction 


Recently there has been a surge of interest in the use of continuum models to study 
the dynamics, stability, and control of vibrating structures. This trend has been 
especially apparent with regard to large flexible spacecraft. Often an important tool 
in these studies is the structure’s natural vibrational frequencies and corresponding 
mode shapes. 

A class of structures which is playing an ever increasing role in these studies is 
the one which consists of flexible beams with tip bodies. We have investigated the use 
of piecewise polynomial spline-based Rayleigh-Ritz-Galerkin schemes for the 
approximation of the natural modes for the transverse vibration of structures of this 


type. 

It is often the case in engineering practice that in using the Rayleigh-Ritz method 
to determine the natural modes of a complex structure the approximating subspaces of 
trial functions are chosen as the span of a finite number of mode shapes for a related, 
but more easily analyzed, structure. In our investigation, which was motivated by one 
such instance of this (see [18]), we observed that the spline-based schemes were 
computationally attractive and yielded significant increases in accuracy and * 
stability. In addition these benefits were achieved for relatively low orders of 
approximation. 

Although we treat a relatively specific structure, it is not difficult to see how 


our results would extend to a more general class of problems. Mathematically, we have 
proven few new results, but rather have modified and applied several existing ones 


which have appeared throughout the literature. The system which describes the 


transverse vibration of a flexible beam with tip body Is a hybrid of ordinary and 


partial differential equations. In section 2 we call upon the well known product space 


theory for functional differential equations to develop an appropriate abstract 


operator theoretic formulation for the eigenvalue problem and to establish relevant 
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spectral results. In section 3, the Rayleigh-Ritz method and the associated 
convergence theory for the approximation of eigenvalues and eigenfunctions of elliptic 
differential operators as developed in [2], [3], and [14] are discussed in the context 
of the present problem. In section 4 we outline the application of these results to 
spline-based methods in the spirit of the treatment in [15]. In section 5 an example 
and numerical results are discussed. 

Notation is standard throughout. The Sobolev space of functions 6 defined on the 

k-1 k 

interval (a,b) for which D <f> is absolutely continuous with D e is denoted by 

H^(a,b). The usual Sobolev inner products and norms are denoted by and 1*!^ 

respectively. The spectrum, point spectrum, and continuous spectrum of a linear 

operator T are denoted by l(T), T, (T), and 1 (T). For A e P(T), the resolvent set 

P ^ 

of T, the resolvent of T at X is denoted by R^(T). 


2. Formulation of the eigenvalue problem and spectral results 

We consider a long, slender, flexible beam of length £ having spatially varying 
linear mass density p and flexural stiffness El which is clamped at one end and free 
at the other with a rigidly attached tip body. The mass properties of the tip body are 
assumed to be: mass m centered at a distance c from the tip of the beam and directed 

along the longitudinal axis of the beam and moment of inertia J about its center of 
mass (see Fig. 2.1). 
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Figure 2 . 1 


Letting u(t,x) denote the vertical displacement of the beam at time t > 0 at 
position x e [0 , £ ] and assuming small deformations (|u(t,x)| « £), the Euler- 
Bernoulli theory for the transverse vibration of a flexible beam yields the equation 
(see [4], [18]) 

(2.1) p(x)D^u(t,x) + D*EI(x)D*u(t,x) - D o(x)D u(t,x) = 0 

where a denotes the internal tension which results from externally applied, temporally 
invariant, axially directed loading. Elementary Newtonian mechanics can be used to 
derive the boundary conditions which describe the motion of the end of the beam with 
the tip body. Translational and rotational equilibrium respectively yield the two 
boundary conditions at x = Z given by (see [4], [18]) 
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(2.2) raD 2 u(t,£) + mcD 2 D u(t,£) - D EI(H)D 2 u(t,)l) + a( Jt)D u(t,£) = 0 

L LX X X X 

and 

(2.3) mcD?u(t,£) + ( J+mc 2 )D 2 D u(t ,£) + ElU)D 2 u(t,*) + caU)D u(t ,J.) = 0. 

L LX X X 

At the clamped end, x = 0, we have zero displacement and zero slope as given by the 
geometric boundary conditions 

(2.4) u(t,0) = 0 

and 

(2.5) D x u(t,0) = 0 

respectively. 

Remark If the center of mass of the tip body were not directed along the longitudinal 
axis of the beam, that is, the tip body had a non-zero mass center offset, the 
resulting set of equations would be nonhomogeneous . Consequently a classical 
eigenvalue problem for the natural modes of vibration would not result. We do note, 
however, that by employing a suitable transformation of the parameters c and J the 
homogeneous part of the resulting equations can be put in the form of (2.1)-(2.5) 

(See [18]). 

The natural frequencies and mode shapes for the system described above are 

determined by assuming a solution to (2. 1)— (2.5) of the form u(t,x) = e itdt <J>(x). The 

2 

eigenvalue problem on <)> with eigenvalues X = w 

D 2 EI(x)D 2 <J,(x) " Da(x)D<f>(x) = Xp(x)<{.(x) 


( 2 . 6 ) 
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(2.7) -DEI( A)D 2 4>(£) + a(OD<f>U) = X(m<(.(jl) + racD<t>U)) 

(2.8) ElU)D 2 <j>(£) + ccr(£)D<}>(£) = X(rac<J>U) + ( J-Hic 2 )D<|>(jO ) 

(2.9) <f>(0) = 0 

(2.10) D<fr(0) - 0 

results* In order to study the eigenvalue problem (2. 6)-(2. 10) , we develop an 

appropriate abstract operator theoretic formulation. 

2 0 

Define the Hilbert space H = R x H (0,£) with inner product 
<(n,<f>), (5.40> H = n T C + 

We make the standing assumptions 


(HI) 

p e C(0,£), p(x) > 0 x e [0, 

a] 

(H2) 

El £ C 2 (0,£), EI(x) > 0 x e 

[0,4] 

(H3) 

a e chOjJl) 



and of course that m,J,c > 0. The smoothness requirements in (HI), (H2), and (H3) 
above can in fact be relaxed. This will be discussed further at the end of the 
section. 

Define the operators M : H ->■ H, Aq s Dom(A) c H -»■ H and B^ : Dom(B) c H + H by 
M(n»«j)) = (MqTi, p<}> ) , 
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m me 

M ° = T+m 2 

L nic J+mc J , 

Dom(A) = { (n , 4>) e H : <(» e H 4 (0,£), <f>(0) = D^> (0) = 0, n = ( 4>( Jt) , D4>(£)) T }> 
A Q ^ = ( (— DEI(£)D 2 <J>(£) , EI(£)D 2 4)(£)) T , D 2 EID 2 <f>) , j = ((<f>U), D<|>(£)) T , <|>), 
Dom(B) = {(n,<f>) e H : <}> e H 2 (0,£), <J>(0) = D<f>(0) = 0, X) = (<f>U), D<f>U)) T }, 
B 0 i = ((a(il)D<(,(£), ca(£)D<f>(£)) T , - DaD<j>). 


Defining Lq : Dom(L) c H H by Dom(L) = Dom(A), = Aq + Bq we consider the 
abstract formulation of (2.6)-(2.10) given by 

(2.11) Lq 4> = XM<J) <{> e Dom(L). 


Since M is invertible, the abstract generalized eigenvalue problem (2.11) is equivalent 
to 


(2.12) L<}) = l e Dom(L) 


where 


L = M Kq = M ^ + M _1 B q = A + B, 


M 1 (n » 4>) = (M Q ri, 1 <J>)> 


with 
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i - 1_ r *- 2 

0 ~ Jm L -me 


-me 

in 




We investigate the spectral properties of the operator L by first characterizing 
the spectrum of A and then treating B as an A-bounded perturbation. Defining 
the <*,»>^ equivalent inner product on H, <«,*>^ by 


it is not difficult to argue that A is densely defined, self adjoint, and positive with 

respect to the <»,»>^ inner product. Consequently A is closed and its 

spectrum, E(A), is real with E(A) c (0,») and E(A) = £ (A) U T. (A). 

P C 

1 2 

For each X e C let <f> and <(> denote two linearly independent solutions to 

A A 

— D 2 EID 2 <}) - A<f> = 0 
P 

<K0) = o 

D<)>(0) = 0 

2tc2 

and define A^ e C by 

-DEI(£)D 2 <f>*(£) — DEI(£)D 2 <f> 2 (£) 

EI(£)D 2 <|>*(£) EI( £)D 2 <J> 2 ( £) 


(2.13) A X = M Q 1 


♦ J(£) ^U) 

D<J>^(£) D<(. 2 (£) J . 


Using arguments similar to those which can be found in [21], it can he demonstrated 
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that A - XI is injective/sur jective if and only if i is injective/surjective. 

A 

Consequently, since A^ is finite dimensional we have the following theorem. 

Theorem 2. 1 The spectrum of A is discrete. That is T,( A) = ?: (A) with X e £(A) if and 
only if det A.. = 0. 

A 

Moreover, A has compact resolvent. 

Theorem 2.2 For each X e P(A) , R^(A) is a compact linear operator. 

Pf 

Let X e P(A) (that is det A^ * 0) and suppose that X is not an eigenvalue of the 
clamped-clamped beam vibration problem 

(2.14) - D 2 EID 2 <(> = X<{. 

P 

(2.15) (f)(0) = D (f)(0) = <t>(Z) = D<t>(Z) = 0. 

Furthermore, let G(x,y; X), 0 < x, y < £ denote the Green’s function corresponding to 
(2.14), (2.15), and X (see [5]). The operator R^(A) can then be written as the sum 
of two bounded operators, 

V A) = T 1 + V 

where T 2 is given by 


T 2 (C,i|0 = ((0,0) T , /qG( • >9 ; X)i|»(0)d0) 


and Tj has range in the finite dimensional space 
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{(n,<J>) e H : n = (<t>00, D<t>U)) T , 4> = + c 2 *\ } - 

Now T 2 is essentially an integral operator on H^(0,£) and as such is compact (see 
[7]). The finite rank operator is compact from which we may conclude that R X (A) is 
compact being the sum of two compact operators. 

Having now shown that R X (A) is compact for at least one X e P(A) , we may 
conclude that A is discrete and hence that (A) is compact for all X e P(A) (see 
[7], Vol III, p. 2291). 

A straightforward application of the spectral theorem for comoact self adjoint 
operators in a Hilbert space (see [20]) yields the following well known result. 

Theorem 2.3 1(A) = {X k ) k _^ with 0 < X^ < X 2 < X^ •••• <+«» and X k +» as 
k -»■ < 0 . Each X^ is of multiplicity 1 or 2 and the corresponding set of 

a m 

eigenvectors, {<}>,}._, c Dom (A), is complete and orthonormal with respect to the 

1C K— A 

inner product. We have 


A A A A A 

<f> = I <$> <t> e Dom (A) 

k=l 

and if X e P(A) 

A 

00 <(n »d>) » <t> k > M -s 

R x (A)(n,<(.) = I 

k=l k 

(tu<(>) e H. 

The operator B is densely defined with Dora (A) c Dora (B) and symmetric (B c B ) 
with respect to the inner product. Furthermore, hypotheses (HI) - (H3) and the 
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Sobolev interpolation inequalities for intermediate derivatives (see [1]) can be used 
to argue that B is A-bounded with A-bound 0 (see [10]). That is, for all 
e > 0 sufficiently small there exist constants a and 6, independent of e, for which 

(2-16) |B+| m < ae _1 UI M + Pe|a{| m » ♦ e Dom (A). 


It follows from (2.16) that L = A + B is self adjoint (see [10]) and therfore closed. 


Since A is self adjoint and positive it is m— accretive and as such satisfies 


and 


I r -x (a) Im 


1 

Re X 


I ar -x (a) Im < 1 


for all X e C with Re X > 0. These estimates imply that 

_i -1 

ae |R -X (A)| M + 8 «|AR. x CA>| m < + 6e < 1 

for e chosen sufficiently small and all X e C with ReX > 0 sufficiently large. We 
may conclude therefore (see [10], page 214, Theorem 3.17) that there exists a 
constant y such that (X e C : ReX < y} c P(L) with R X (L), X e P(L), compact. 

It is not difficult to argue that 


(2.17) y > -4aB. 

As was the case with A, the spectral properties of L are easily characterized 
using the spectral theorem for compact self adjoint operators on a Hilbert space. 
Indeed, the conclusions of Theorem 2.3 are valid with A replaced by L and with the 
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exception of the fact that the spectrum of L is bounded below by -4ct6 rather 
than 0. 

The smoothness assumptions in (HI) - (H3) can be relaxed to p, El, a e 1^(0, fc) 
with El bounded away from 0 and a defined at Z by turning to a weak formulation. 
Define the space 

V = { ( n , d>) e H : (j> e H 2 (0,£), <K0) = D<f> (0) = 0, ti = (<f>U)> b<|>U)) T } 
with inner product 

<i r i 2 > v - < 1)2< f , 1 * ° 2 <f> 2 > o 

where ^ = ( ( <t> ± ( «.) , D^U)) 1 . (f^), i = 1,2. The usual dense embeddings Vc H e V' 
hold with the injection V <*H compact. Define the bilinear forms on V x V 

a(i, i) = <EID 2 <t>, D 2 i|>> 0 

ca(£)Dc(>(£)Di|;(£) + <aD<f) , DiJ)>q 

A /\ A A 

a(<{>,4>) + b(<j> ,i|>) . 

The form a(*,») is V-elliptic (see [17]) while the form 1(»,*) is V-H-elliptic. 
Indeed recalling (2.17) we have 

!(♦.♦) + ^ag|}| 2 > 6|*|J 

for $ e V. The forms a(»,»), b(».»), and !(•,•) are related to the operators A, B, 


b(<t>,4») = 

and 

i(J,i) = 
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and L defined earlier in the usual manner. Consider the form a(«,»). Using the Riesz 
Theorem define the operator A^ : Dom (A^) c h -v H by 

A A A A 

Dom (A^ ) = {$ e V : The map iJj -»• a(<f>,iji) is continuous on 

V with respect to the | • | norm on H> , 


a( <}>>'!') = <f> e Dora (A^), ^ e V. 

It can be shown that 

Dom (Ap = {} e V : EID 2 <f> e H 2 (0,£)> 

and if El is sufficiently smooth, that Dom (Aj) = Dom (A) and Aj = A. In addition if 
the inner product is interpreted as the duality pairing between V and V' then 

A 1 can be extended to an operator A^ e L(V,V'). Similar correspondences exist between 
the form b(*,») and the operator B and the form 1(.,.) and the operator L. 

The weak form of the eigenvalue problem (2.12) is given by 

(2.18) 1(<}>,iJj) = ip e V. 

Standard results (see [17], page 78, Corollary 7D) yield the existence of a set of 
orthonormal (with respect to the inner product) eigenvectors 

tx> 

^k^k=l c Dom ^ A l^ c V which are complete in H and a corresponding set of eigenvalues 
1 which satisfy —4af5 K A^ < A 2 •••• < «> with A^ -*■ -H» as k -»• <». 
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3. Approximation and Convergence - The Rayleigh-Ritz - Galerkin Method 

Since the addition of a constant multiple of the identity to an operator simply 
translates the spectrum by the same constant value, we assume throughout, without loss 
of generality so far as approximation is concerned, that L is positive. For (2.17) 
implies that the operator = L + 4aBI satisfies 


<L X J,i> M = <04 ,J> m + 4aB<i,i> M > -4aB<J,^> M + 4aB<£,$> M = 0. 


The Rayleigh-Ritz method for the approximation of the eigenvalues and eigenvectors of a 
positive self adjoint operator is based upon their characterization as the extrema and 

A, A 

critical values (stationary points) of the Rayleigh quotient. For 6 e V, <J> * 0 define 


Rui - Ki, ♦)/<;, ;> M . 


^ A ^ A A A 

We note that for <)> e Dom (L) and t|> e V we have l(<f>,^) = and that V is the 

A A A 2 

closure of Dom(L) in V with respect to norm = (1 (<!>j<I>)) on V* The 

00 00 

eigenvalues and corresponding eigenvectors can then be characterized 

by (see [9], [22]) 

(3.1) = min{R[<()] : <j> e V, <|> ^ 0, <<}> »<(>..)> ^ = 0» j = l,2,...,k— 1} — R[<t>^]. 

There exist other characterizations of the spectrum equivalent to (3.1) which are often 
more useful for computational and theoretical purposes (see [14]). 
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For each N - 1,2,... let V N denote a finite dimensional subspace of H with 

N 

V N c V, dim V N = k N and V = span The classical Rayleigh-Ritz method 

consists of determining the extrema and critical values of R[.] over V N . Requiring 
that the first variation of R[»] vanish over V N results in the matrix generalized 
eigenvalue problem 

(3.2) L N v N = xW 


where 


FT N i 1/?N 7N. 

[L ] ± j = l(i p ±9 


i>j — 1,2, ... ,k 


N 


and 


[M N ] = <JJ, i,j = 1,2,.., k N . 


The matrices L N and M N are real, symmetric, and positive definite. Let 


C N = (H N ) _ 2 L N (m N } " 2 # 


Then the generalized eigenvalue problem (3.2) is equivalent to the standard eigenvalue 
problem 


JN N N 
C w = X w 


JL_ 

where w N = (M N ) 2 v N . The matrix C N is real, symmetric, and positive definite. Tt has 

N 

k N positive eigenvalues U^}£ =1 > 0 < x” < x!J < ••• < X N N and k N corresponding 

k 
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N,k‘ 


N 


eigenvectors which can be chosen to he orthonormal, 


, N n T N . , . , 0 v N 

(w i ) Wj = 6^ i,J = l,2,...,k . 

N k N *N k N 

The are known as Reyleigh-Ritz approximate eigenvalues and , 

k N 

r N A N 

4> . = l (v. ),ij>, the corresponding Rayleigh-Ritz approximate eigenvectors of L. The 
J i=l J 1 


?N,k 


N 


^k^k=l are ort ^ onornia l H with respect to the <*,*>^ inner product, 


,^N *N n , N,T m N N / NvT N _ . 

<4>i > <frj> M = (v ± ) M Vj = (w^ Wj = 6 ± j 

i,j = l,2,...k N . 

The weak formulation (2.18) is the basis for the Galerkin method. It leads to the 

A N N N 

finite dimensional eigenvalue problem which consists of finding <f> e V and X that 
satisfy 


o v w a n 7N, . n a n ;n n ;n ^ u n 

(3.3) 1(<|> , ) = X <4> , i|> > M , e V . 


For the problem considered here, the Galerkin equations (3.3) and the Rayleigh quotient 
characterization (3.1) both lead to the same matrix generalized eigenvalue problem. 

The convergence theory for the Rayleigh-Ritz method for the approximation of the 
spectrum of self adjoint elliptic differential operators is well documented in the 
literature. Convergence results for the eigenvalues and eigenvectors of Sturm- 
Liouville systems can be found in [2]. These ideas were then extended to apply to more 
general problems and improved upon in [3] and [14]. Although the operators Lq and M 
are not precisely of the form of those considered, it is not difficult, however, to 
argue that the formulation we have employed renders the existing convergence theory 


directly applicable 
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N N m 

Let P : V V denote the orthogonal projection of V onto V w with respect to 

A A A A 

the <*>*> L inner product on V defined by = l((f>,ij>)« Assume that the subspaces 

V N have the property that 


lim I P N i - i\ = 0 i e V. 

N-»-« 


Theorem 3.1 (Eigenvalue Error Estimate). Let k be a fixed positive integer with 
k < k . Then the sequence {X^} converges to from above, there exists a 
positive number N such that for all N > N 


\ \ + 1 




a-r ; hVi'S* 2 

J=1 


and hence 


< K + v 


I \v\ 

j-1 3 



with 



as 


N 4- «• 


Theorem 3»2 (Eigenvector Error Estimates), Suppose X fc is of multiplicity 

N 

m^ + 1, m^ = 0 or 1 and k + m^ < k • Then: 

(i) If m^ = 0 (X^. is simple) there exist positive numbers y and N depending only uoon 
k for which 

U k - * t l L < y|p - + k l L 
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for all N > N. 

/v 

(ii) If in^ = 1 (X k is of multiplicity 2), let <f> be any element in 

A A IT ^ 

span{£, , $ } and let $ be the orthogonal projection of <f> onto 

tC rCr 1 

span{iJJ , with respect to the innerproduct. Then, there exist 

positive numbers y and N depending only upon k for which 

— k+1 

|J N - Jl L < yUIm l 

j-k 

for all N > N. 


In the next section we describe a particular class of schemes which are based upon 
choosing the V N as subspaces of spline functions. We use Theorems 3.1 and 3.2 together 
with the approximation properties of splines to establish convergence and estimate 
rates of convergence* 


4 Spline Approximations 

N N 

For each N = 1,2 ... let tt denote a partition of [0,£], tt : 0 = x Q < xj < x 2 

... < x M = £ and let h N = max (x - - x^) and g N = min (x ± - x^). 

N i=l , 2 , . . . ,N 1 i=l,..,n 

For ra = 1,2 ... let S(ra, tt, z) denote the order 2m (degree 2m - 1) polynomial spline 

N i ,jN— 1 

spaces corresponding to the partition tt and the node incidence vector z e 

(see [19]). Recall that S(m,u N , z) c H 2m y (0,£) where 1 < V = raax z i < m 

N-l i=l,2,...,N-l 

and dim S(m, tt N , z) = 2m + n(N, z) where n(N, z) = V z. . 

i=l 

Of particular interest in practice are the following two special cases: 



- 18 - 


(1) z (m,m,...,m). In this case S(m,ir , z) - H^(it ), the Herroite spaces, 

dim H ( tt N ) = ra(N+l) and H (tt N ) c H m (0,£). 
in m 

(2) z = (1,1, . .. ,1). In this case S(m,7r N ,z) = S^ir 1 *), the standard polynomial spline 
spaces, dim S m (r N ) = N - 1 + 2ra and (tt N ) c H 2m-1 (0,£). 

For m > 2 let S (m,ir ,z) = {s e S(ra,ir ,z); s(0) = Ds(0) = 0} and define 

V m ( z ) = {<}> = ((<f>(£), (£ ) ) "^ , ) : 4> e S^(m,ir^, z)} . 

Then k (z) = dim V (z) = 2(m-l) + n(N,z) and since m > 2, V^(z) c V. For k > 2 
ui m rn 

define 

V k = { (ti , C f>) e H: <J> e H k (0,£), (f> (0) = Dd>(0) = 0, n = (MO, D<f>(0) T }. 

Then V k c V, V k+1 c V fc , k = 2,3,... and V 2 = V. 

In order to simplify notation, we suppress showing explicit dependence on z with 

the understanding that the spaces and the interpolation and projection operators 

defined below do in fact depend upon the choice of the node incidence vector. 

Define the interpolation operators I N , : H k (0,£) •+• S(m,ir N , z) by letting 

in j 1 

N N 

I m,k < * > denote the unique element in S(ra,r z) (see [19]) which satisfies 

D ^(<J> “ 1^ k <f>)( x ^) = 0 0 < 1 < min(k-l, z^-1) 

0 < i < N 

I^k^^i) = 0 min(k-l , Zj-1) < j < z i ~l 

0 < i < N 

where z n = z„ = ra. For ra > 2 define , : V, -»• V N by 

D N m,k k m 3 
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m k* = «♦![(*>» V!U)) T , A 

m,k T mm m 


where <|> = ((<}>(£.) , 0<t>(O) T , <t>) e and <f>^ - 

N °° 

If the sequence of partitions {tt satisfies the uniformity condition 

h^ ^ 

— < x for some x > 1 independent of N and lim h = 0 we have the following error 

g N N+« 

estimates. 


A T 

Theorem 4» 1 For m = 2,3,.,,, k = 2,3, 2m and (f> = (($(£), e 

(i) " il V < TCH N ) k - 2 |D k 4>| 0 

and 

(ii) 2 4*' “ <f*l v ♦ 0 as N °°> l e v 2 = V 
where y is a constant independent of N, 


Pf 

Statement (i) is a direct consequence of standard spline interpolation error 

A 

estimates (see [19]). To argue the validity of (ii) let <j> e V. Then for 
\|> e H^(0 , i) 


S'* J - 

m, 2 V 


l v < 1^ - *l 2 < I *1 - *Jl 2 + " ^2 + I* “ *•: 


where = 1^ 0 iJj. Now (see [19]) 
m m, L 


^m "O 2 < Y ' D ^ (I m,2 (<!> “ ♦)>!<> < - ^'0 < " ^2 


and therefore 
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(4.1) I 2^ “ ^ I v < ^ + Y)|<J> ~ ’M 2 + I’J'm “ ,f ’l2* 

Corresponding to z and 1^ . = 1^ (z) there exists a z and I N _ = I N „(z) such that 

th, J to, j 

I m,2^ = ^,3*- Therefore 

(4.2) |/ - *| 2 = |I^ 2 * - if) | 2 = - *| 2 < Yl h N |D\| 0 . 

The estimates (4.1) and (4.2) together with the density of H 3 (0,O in H 2 (0,O and the 
N 

fact that lim h = 0 yield the desired result. 

N-Voo 

Let P m : V "*■ V m denote the orthogonal projection of V onto computed with 
respect to the < # »’> L inner product. Arguments similar to those found in [15] 
together with the interpolation error estimates given in Theorem 4. 1 can be used to 
verify 

Theorem 4.2 If l = ((<j>U), D<J>U)) T ,<f>) e V k with 2 < k < 2ra then 

(i) l P 5 " *I M < Y(h N ) k |D k <H Q , 

^d) |P^ <J> - (J>|^ < y(h^) 2 |D k <J>| 0 and |P^J <() - <t>|^ + 0 as N + » if k=2, 

(iii) |D J (^ - *)| 0 < Y (h N ) k " j |D%| 0> j-0, 1, 2 and 

l° 2 (<f>^ ■ OIq + 0 as N + » if k=2, 
k-j- — 

(iv) |D j (^ - ^ < Y (h N ) 2 |D%| 0 , j-0,1, 

N ^ N N T N 

where p m <f> = ((<f> m (£)> Dt|> m (£)) , <}) jn ) and y is a constant which is independent of 
<(> and N. 
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Theorems 3.1, 3.2 and 4.2 yield the following eigenvalue and eigenvector error 
estimates (see [15]). 

a tti A N m 

Theorem 4.3 Suppose X ^ is simple and <f>^ e V^. Let X^* and denote respectively 

rh N 

the i tn Rayleigh-Ritz approximate eigenvalue and eigenvector for L computed over 

N N 

where N and ra have been chosen so that k = dim V > i. Then 

m m 

(i) 0 < X N ’ m - A ± < y(h N ) 2 ^ r ” 2 ^ and - X ± + 0 + as N + « if r=2, 

(ii) " 4> ± | L < y(h N ) r 2 |^ r <t> i l 0 and + 0 as N + » if r=2, 

(iii) i^ ra -i i i H <Y(hViD r ^ i i 0 , 

(iv) |D J (*J ,m - cf. ± ) | Q < Y(h N ) r_1 |D r <t 1 l 0 , .1=0, 1,2 and 

| D 2 (<J>^’ m - + 0 as N + » if r=2, 

r-j- I 

(v) | - <f> i ) I „ < y(h N ) 2 |D r <(> i l 0> j =0>1 

where i ± = ((4» ± (jO, D^U)) 1 , * ± ), JJ*® - D<|,J’ m U)) T , <f»^’ m ) , y is a 

constant which is independent of N and r=rain(k, 2m). If X^ is of multiplicity 2, we 

A A A l^[ ftl A iq" ^ A 

have (i) - (v) above holding with <j>^ replaced by <|> and 6^* by <f> * where $ is any 

A A A N m A 

element in spanf^, and $ 9 is the orthogonal projection of <fr onto 

A N m m 

span{<^* , computed with respect to the < • , • inner product. 

N N N 

Remark When V is defined using S (tt ) with tt a uniform partition the L error 
m m co 

estimates given in the previous two theorems can be improved. Orders of convergence in 
Theorem 4 # 2 (iv) and 4.3 (v) can be increased by ^ • 
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5. An Example 


We consider the approximation of the natural frequencies for the transverse 
vibration of a uniform, cantilevered beam with tip body subiect to an axially directed 
base thrust. More specifically, we assume that the linear mass density p and flexural 
stiffness El are spatially invariant. Furthermore we assume that the beam is acted 
upon by a constant, axially directed load Pq > 0 applied to its base inducing an 
acceleration in the positive x-direction (see fig 2.1). It can be derived (see [181) 
that the internal tension a is given by 


0(x > - - p o( 


m + / Q p 


m+£p 


We note that the formulation described above leads naturally to a classical 

buckling problem. Indeed, for Pq sufficiently large the beam will buckle. The critical 
cr 

buckling loads, Pq are those values of Pq for which steady state solutions to (2.1) - 
(2.5) exist; that is those values of Pq for which X = 0 is a solution to the 
eigenvalue problem (2.6) - (2.10). For the problem stated here, the critical buckling 
loads can be computed directly as the roots of a transcendental equation involving 
Bessel functions of the first kind (see [18]). 

The B-spline representation for polynomial spline spaces leads to computationally 
efficient algorithms for the evaluation of the spline functions, their derivatives, and 

inner products (see [16]). We considered a scheme using the standard polynomial spline 

Nv N N 

spaces S^tt ) defined over the uniform partition tt with h" = £/N. In this case the 

relevant formulas become especially simple. 

For k > 0 let (x-y)^ = (x-y) k (x-y)^ where 
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(x-y) 


0 

+ 



x < y 
x > y . 


Define the fundamental B-spline of order 2m (degree 2m-l) with knots {0,1,2... 2m} and 
support [0, 2ra] by 

„ f N .2m, N 2m-1, , lN i r 2mN / ,v2m-l 

V x) ’ S y <x - y) + ly-0 ' * <-» r 

where 6 is the forward difference operator on y; 6yf(x,y) = f(x,v+l) - f(x,v). If 
we let 


B^ ,m (x) = (x - (i-m)h 1 *)), i=l-ra, . . . ,N+m-l 

then S (ir^) = span{B^’ m }^ + ^ 1 ^ The basis splines B^ ,ri for S (ir^) are scaled-- 
m i i=l-m. i ra 

translates of the fundamental B-spline, B m , and as such have the desirable property 
that their supports, ((i-m)h^, (i-hn)h^) are local. A. basis for 

S^(tt^) = {s e S (tt N ) : s(0) = Ds(0) = 0} is easily constructed from the B' I,ra by 
mm 1 

taking linear combinations of basis elements which do not satisfy the boundary 
conditions at zero to form basis elements which do. We obtain 

S^(ir^) = s P an ^i * ^ or exara P^- e » t ^ ie case m= 2 one nossible choice for 

the is given by 


6 


N * 2 = Bq’ 2 - 2B^’ 2 - 2B^’ 2 


e 


N ’ 2 = Bj’ 2 , j = 2,3, .. . ,N+1. 



-24- 


Defining = ((Bj ,ni (£), DB^’ ra (£)) T , 6^ ,n ) e V, j = l-m+2 , . . . ,N+ra-l , m > 2, we 

have = span {8^* ^j=l-m+2 * matr i ces and which appear in the generalized 

eigenvalue problem (3.2) are given by 


(5.1) [h\^ = El /J D 2 B^ m D 2 B j N,m - ^ /J (mf(£-x)p)DB^ m DB^ ra 

ZZT- ^ ,m coDeJ* TO (i) 


and 





— — 

(5.2) 

ij - <’"(*)! [” c 

me j 
J+tncl 

e"' m U) 

i»j = 

1 -m-f 2 ) • • • ) N+m- 1 • 




£ Q N >m 0 N,m 


+ P W ^ 


Since the B-splines are polynomials of degree 2m- 1 on each subinterval of tt N , the 
integrals in the above expressions can be computed exactly using a composite Gaussian 
Quadrature formula. The efficient evaluation of the B-splines and their derivatives is 
facilitated by the fact that the fundamental B— spline satisfies the recurrence 
relations 


B m (x) = x 2 B m _ l( x) + (2x(2ra-x) - 2m)B m _ 1 (x-l) + (2ra-x) Z B (x-2) 
DB m^ x ^ = (2m-l) [xB m _j(x) + 2(ra-x)B m _^ (x-1) - (2m-x) B m _^ (x-2) ] 
D 2 B m (x) = (2m-l)(2m-2)[B m _ 1 (x) - 211^ (x-1) + B^^x-2)]. 
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The local support property of the B-splines necessarily implies that the matrices 
and M N are banded. This can be exploited when (3.2) is solved (see [6]). 

We chose £ = 1*0, p = 2,0, El = 1*0, m = 4*0, c = *1, J = .52 and Pq = 2*5 and 

solved the generalized eigenvalue problem (3.2) with the matrices and given by 

(5.1) and (5.2) respectively for ra = 2 (cubic splines), 3 (quintic splines) and 4 

(septic splines) and various values of N. The eigenvalue problem (3.2) was solved 

using the IMSL routine EIGZS. This routine finds the Cholesky decomposition R N (R N ) T of 

M N and then uses the QR algorithm to compute the eigenvalues of the transformed matrix 

(R^) - ^L^(R^) -r *‘ . Our results are given in Tables 5.2, 5.3, and 5.4. Recall 

that dim(S (tt N )) = N+2ra-3. 
m 

We also computed approximate eigenvalues using the Rayleigh—Ritz method over two 
modal subspaces. The spaces were taken to be the span of the first N natural mode 
shapes for the corresponding (1) claraped-f ree beam (m = c = J = Pq = 0) and (2) 
clamped-free beam with tip body (Pq = 0). These results are given in Tables 5.5 and 
5.6 respectively. 

As an independent basis for comparison, we also computed approximate eigenvalues 
by using an iterative root finding technique to compute the zeros of the frequency 
equation corresponding to (2.6) - (2.10). A numerical integrator was used to find two 
appropriate linearly independent solutions to the initial value problem (2.6), (2.9), 
(2.10) for a given value of X in each iteration. The frequency equation arises from 
the setting to zero of the determinant of a matrix which is similar in form to given 
in (2.13) (see [18]). The first five eigenvalues computed in this manner are given in 


Table 5. 1 below. 
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i 

X i 

i 

.0365798 

2 

9.597101 

3 

271.624402 

4 

1934.262208 

5 

7352.906500 


Table 5.1 


U 

N+ 2 

4 

8 

16 

32 

1 

.0367460 

.0365892 

.0365804 

.0365799 

.0365799 

2 

9.611071 

9.597834 

9.597144 

9.597104 

9.597107 

3 

279.748161 

272.659779 

271.677636 

271.627567 

271.624598 

4 


2008.349824 

1937.362657 

1934.430532 

1934.272355 

5 


7808.026186 

7406.496727 

7355.451937 

7353.054424 



Table 

5.2 - Spline Scheme - m = 2 



+i 

N+ 2 

4 

8 

16 

32 

1 

.0365798 

.0365798 

.0385798 

.0365798 

.0365798 

2 

9.597110 

9.597101 

9.597101 

9.597101 

9.597101 

3 

271.637063 

271.626382 

271.624408 

271.624403 

271.624403 

4 

2024.194063 

1937.055444 

1934.264769 

1934.262236 

1934.262227 

5 

7730.251147 

7358.613802 

7353.153706 

7352.907163 

7352.906727 


Table 5.3 — Spline Scheme — ra = 3 
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u 

N+ 2 

4 

8 

16 

1 

.0365798 

.0365798 

.0365798 

.0365798 

2 

9.597101 

9.597101 

9.597101 

9.597101 

3 

271.624474 

271.624406 

271.624403 

271.624392 

4 

1935.739405 

1934.300633 

1934.262216 

1934.262210 

5 

7362.669220 

7353.063985 

7352.907860 

7352.906573 



Table 5.4 - Spline Scheme - m - 

4 


4-i N> 

2 

4 

8 

16 

1 

.0438998 

.0391454 

.0375638 

.0370661 

2 

20.580764 

13.553588 

11.326133 

10.581254 

3 


445.668737 

340.718106 

310.220838 

4 


3401.113918 

2441.500175 

2211.540908 

5 



9373.564654 

8425.737638 

Table 5. 5 - Clamped Free Modes 


U N+ 

2 

4 

8 

16 


i 

.0382750 

.0365966 

.0365801 

.0365608 


2 

9.609008 

9.597279 

9.597107 

9.571124 


3 


271.642143 

271.624818 

270.693532 


4 


1934.288529 

1934.263452 

1927.899003 


5 



7352.912644 

7328.632345 



Table 5.6 - Clamped Free with Tip Body Modes 
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Frora the tables above it is immediately clear that the spline schemes exhibit 

rapid convergence with no apparent deterioration of accuracy for large N. Rapid 

Convergence was also observed in higher frequencies (i > 5) than those shown in the 

tables. The cantilever mode based scheme converged very slowly. The cantilever with 

tip body modes performed better. However, there is an apparent stability problem when 

N becomes large. This is most likely a consequence of the high frequencies involved 

(which must be determined as the roots of the transcendental equation det A, = 0) and 

A 

the fact that the stiffness matrices L N are full. The cantilever with tip body modes 

performed optimally when N - 12. Rapid eigenvalue convergence was also observed with 

Hermite spline-based schemes. However these schemes did exhibit some stability 

problems when N was large. Their overall performance was inferior to 
N 

the S (tt ) - based methods, 
m 

Finally we note that if damping is introduced either into the beam or tip body 
dynamics it is likely that the resulting stiffness operator would not be self 
adjoint. Consequently the theory presented here would no longer be directly 
applicable. However, the results in [8] and [13] or [11] and [12] may well apply 
depending upon whether or not the associated resolvent operators were compact. We have 
not as yet considered the damped problem. 
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